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Abstract —  Although  magnetic  resonance  (MR)  tagging 
has  been  shown  to  be  a  useful  tool  in  myocardial  mo¬ 
tion  quantification,  its  clinical  utilization  is  limited  as  cur¬ 
rent  available  methods  generally  either  lack  computational 
speed  or  require  extensive  user  intervention.  Recently,  the 
harmonic  phase  imaging  (HARP)  technique  has  been  pro¬ 
posed  to  look  at  the  phase  information  of  the  tagged  images 
[1],  [2].  HARP  imaging  promises  to  overcome  the  limita¬ 
tions  of  existing  methods  in  terms  of  both  computational 
speed  and  automation.  Motivated  by  this  work,  we  present 
mathematical  analysis  providing  a  signal  processing  per¬ 
spective  on  the  HARP  technique.  This  new  perspective 
provides  a  clearer  understanding  of  how  tags  can  be  accu¬ 
rately  tracked  using  highly-filtered  data. 
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I.  Introduction 

MR  tagging  is  a  powerful,  non-invasive  method  to  quan¬ 
tify  myocardial  wall  motion  and  identify  abnormalities  in¬ 
dicative  of  regional  ischemia.  Tags  are  temporary  markers 
created  by  spatially  encoded  saturation  planes  whose  de¬ 
formations  follow  the  motion  in  the  underlying  tissue  [3], 
[4].  After  administration  of  tags  with  spatially  selective 
rf  pulses,  the  deformation  is  observed  through  a  sequence 
of  images  acquired  over  the  cardiac  cycle.  Fig.l.(a)  is  a 
planar  tagged  image  acquired  at  end-diastole  right  after 
the  administration  of  tags  and  Fig.l.(b)  is  acquired  near 
end-systole  showing  the  deformation  of  tags.  Both  images 
have  been  cropped  to  display  a  region  of  interest  (ROI) 
around  the  heart  for  better  visualization  purposes. 

Although  creating  tagged  images  is  a  simple  process, 
extracting  the  tag  deformation  within  the  myocardium  is 
a  serious  challenge,  especially  considering  the  fact  that 
tags  are  only  couple  pixels  wide.  At  each  image,  one  has 
to  trace  myocardial  contours  to  detect  the  end  points  as 
well  as  the  deformed  tags  within  the  myocardium.  The 
complex  motion  of  the  heart  makes  it  difficult  to  utilize  a 
parametric  characterization  of  the  tag  deformation.  Be¬ 
cause  of  short  scan  times,  tagged  images  are  often  poor 
quality,  suffering  from  low  signal-to-noise  ratio  (SNR). 
Furthermore,  T1  recovery  of  the  tagged  tissue  eliminates 
the  use  of  simple  intensity  based  methods.  Despite  all 
these  problems,  there  are  numerous  successful  approaches 
in  the  literature.  Tag  profile  fitting  [5]  and  correlation 
methods  [6]  have  been  successfully  implemented.  Unfor¬ 
tunately,  such  methods  are  computationally  intensive,  re¬ 
quiring  long  processing  times. 

There  are  several  approaches  that  look  into  the  Fourier 
content  of  the  tagged  images  to  address  the  long  compu¬ 
tation  time  problem  [1],  [7].  Especially,  the  HARP  imag- 
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(a)  (b) 

Fig.  1.  Planar  tagged  image  at  (a)  end-diastole,  and  (b)  end-systole 


ing  technique  appears  to  be  very  promising  as  a  clinical 
tool  in  terms  of  speed  and  minimal  user  intervention  [1], 
[2].  In  HARP  imaging,  spatial  modulation  of  magnetiza¬ 
tion  (SPAMM)  tagging  is  defined  as  a  modulation  process. 
Hence,  the  discrete  Fourier  Transform  (DFT)  of  a  tagged 
image  has  spectral  peaks  at  integer  multiples  of  the  tag¬ 
ging  frequency.  Bandpass  filtering  the  proper  peak  results 
in  an  image  in  the  spatial  domain  whose  phase  values  give 
the  tag  deformation. 

In  this  paper,  we  present  a  mathematical  analysis  of  the 
tag  deformation  from  a  signal  processing  perspective,  and 
validate  the  results  of  the  HARP  imaging.  Our  main  con¬ 
tribution  to  the  literature  is  the  introduction  of  tagging 
as  a  sampling  process.  We  believe  that  analyzing  tagging 
from  a  signal  processing  perspective  will  initiate  the  in¬ 
troduction  of  new  methods  to  the  tag  analysis  problem. 
Furthermore,  a  mathematical  derivation  of  the  phase  val¬ 
ues  for  the  tag  points  is  given  in  this  work. 

II.  Methodology 

From  a  signal  processing  perspective,  tagging  can  be 
defined  as  a  sampling  process  such  that  tagged  image  is 
simply  the  difference  of  the  original  image  and  its  sampled 
version  at  the  tag  locations.  Hence,  the  tagged  image  can 
be  expressed  as: 

o(x,  y ,  t )  =  u(x,  y,  t)  -  u(x,  y,  t)tag( x,  y,  t)  (1) 

o(x,  y ,  t )  :  tagged  image 

where  u(x ,  y,  t)  :  untagged  image 

tag(x,y,t )  :  tag  pattern 

In  signal  processing,  a  train  of  impulse  function  lines  is 
called  the  shah  function  ( III(x,y )).  In  an  N  x  N  image, 
the  horizontally  placed  planar  tag  pattern  with  a  tag  sepa¬ 
ration  of  T(t)  can  be  represented  in  terms  of  shah  function 
as  follows: 
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Lt%)J 

tag(x,y,t)  =  Ill(-jj-r)  =  ^  6{y-nT(t))  (2) 

^  '  n= 0 


t  is  the  time,  a;  and  y  are  the  spatial  dimensions  satisfying: 


0  <  x  <  N  —  1  0<y  <N  -1 


Tag  separation,  T(f),  is  defined  as  a  function  of  time  be¬ 
cause  its  value  will  change  for  the  non-stationary  parts  of 
the  tag  pattern  (i.e.  within  myocardium).  The  two  di¬ 
mensional  (2D)  DFT  of  Eq.l  in  spatial  domain  gives  [8]: 

0(u,  v,t)  =  U (it,  v,  t)  —  —U(u,  v,  t)  <8 >jv  Tag(u,  v,  t)  (3) 

where  0  <  it  <  N  —  1  and  0  <  v  <  N  —  1.  ®jv  is  the 
circular  convolution  operator  with  period  N.  The  DFT 
of  a  shah  function  is  an  impulse  train  ( comb  function)  [9]. 
Hence,  Tag(u,v,t)  is  given  by: 

Tag(u,  v,  t)  =  T(t)III  S(u) 

Lt^J  /  N  \ 

=  TW  IZ  6\U,V~kT(f))  ^ 

From  Eq.4  and  Eq.3,  0(it,  v,  t)  then  becomes: 

_jx_ 1 

T(f)  /  N  \ 

0(u,v,t)  =  U(u,v,t)~—  u  \u’v~kJY)J  ^ 

Eq.5  reveals  an  interesting  property  of  tagged  images: 

Property  1:  The  Fourier  spectrum  of  the  tagged  image 
has  replications  of  the  untagged  image ’s  Fourier  spectrum 
separated  by  pixels  from  each  other. 

It  is  important  to  note  that  although  the  Fourier  spec¬ 
trum  of  the  tag  pattern  is  a  comb  function  at  the  time 
of  administration,  it’s  not  exactly  true  for  the  rest  of  the 
frames  because  of  the  tag  deformation  in  time.  Fortu¬ 
nately,  tags  do  not  deform  significantly  from  their  initial 
position,  thus  their  local  frequency  won’t  change  much. 
Furthermore,  tag  pattern  is  applied  to  the  whole  imaging 
volume.  Stationary  tags  dominates  the  pattern  and  forces 
it  to  closely  follow  the  behavior  of  a  shah  function. 

Fig. 2. (a)  is  the  natural  logarithm  of  the  Fourier  spec¬ 
trum  of  the  image  in  Fig.l.(b)  in  which  T( 0)  =  10  pixels. 
The  DC  component  is  shifted  to  the  center  of  the  image 
and  natural  logarithm  is  taken  to  compress  the  dynamic 
range  for  display  purposes.  The  goal  is  to  find  tag(x,y,t) 
and  extract  the  tag  deformation.  Therefore  we  are  inter¬ 
ested  in  finding  T(t).  Bandpass  filtering  the  first  spectral 
peak  as  shown  in  Fig. 2. (a)  gives: 


Ft(u,v,t)  = 


0(u,v,t)  (u,v)  <  passband 

0  elsewhere 


(6) 


(c)  (d) 


Fig.  2.  (a)  Fourier  spectrum  of  Fig.l.(a)  and  the  band-pass  filter 

on  the  first  spectral  peak,  (b)  Unwrapped  phase  magnitude  of  the 
IDFT  of  the  first  spectral  peak  in  (a),  (c)  magnitude  of  the  IDFT 
of  the  first  spectral  peak  in  (a),  (d)  Detected  tag  locations 


The  inverse  DFT  of  H(u ,  v,  t)  is  a  complex  image.  Apply¬ 
ing  a  circular  bandpass  filter  as  in  Fig. 2.  (a)  corresponds 
to  convolving  the  image  with  a  jinc  function  in  spatial 
domain.  Hence,  the  jinc  function  acts  as  an  averaging  fil¬ 
ter  on  the  image  and  creates  a  significant  blurring  in  the 
magnitude  image  as  shown  in  Fig.2.(c).  h(x,y,t)  can  be 
expressed  as: 


T(t)  r  j 2 7r ( iv )  ■ 

h(x,y,t )  =  — ^-jinc  (r(T(t))  Ojv  \u(x,y,t)e  "no  v 


with  its  phase  given  by: 


27 r 


£(h(x,y,t)  =  y  +  ir 


(7) 

(8) 


The  7 r  in  Eq.(8)  comes  from  the  minus  sign,  as  —  1  =  eJ7r. 
One  would  expect  to  have  the  tag  locations  close  to  integer 
multiples  of  tagging  period  T(t).  Hence,  replacing  the  y 
with  kT(t. ),  in  Eq.  8,  the  argument  of  the  phase  values  of 
tag  patterns  can  be  approximated  as: 

arg (l(h(x,y,t))  «  ^L-kT{t)  +  7iy)  mod  (2tt) 

~  7T  (9) 


The  argument  of  the  phase  is  defined  between  [— 7r,7r). 
Because  of  the  tag  deformation  and  the  convolution  with 
the  jinc  function,  the  argument  of  the  phase  values  at  tag 
locations  might  change  sign  and  jump  between  —  7r  and 
7 r.  The  absolute  value  is  taken  to  solve  this  problem.  Be¬ 
cause  the  spectral  peaks  on  each  side  of  the  DC  component 
are  anti-symmetric,  only  one  of  the  first  spectral  peaks  is 
band-pass  filtered  to  eliminate  phase  cancellation.  Either 
of  the  first  spectral  peaks  can  be  band-pass  filtered  as  the 
magnitude  of  the  phase  information  is  used.  Detected  tags 
are  displayed  together  with  the  tagged  image  in  Fig.2.(d). 
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III.  Results 
A.  Computer  Simulations 

A  synthetic  image  with  a  tag  separation  of  ten  pixels  is 
created.  Fig. 3. (a)  shows  the  computer  generated  phantom 
prior  to  deformation.  Sinusoidal  deformations  (AT)  are 
applied  to  the  image  to  obtain  Fig.3.(b): 

AT  =  — ^ sin{4> )  0  <  <j>  <  7r  (10) 


(g)  (h) 


Fig.  3.  (a)  Non-deformed  tag  image,  (b)  Deformed  tag  image,  (c) 

Fourier  spectrum  of  the  non-deformed  tag  image,  (d)  Fourier  spec¬ 
trum  of  the  deformed  tag  image,  (e)  Unwrapped  phase  magnitude 
of  the  IDFT  of  the  first  spectral  peak  in  (d).  (f)  magnitude  of  the 
IDFT  of  the  first  spectral  peak  in  (d).  (g)  Corresponding  tag  points 
detected  from  the  phase  image,  (h)  Tagged  image  and  the  detected 
tags  added  together. 


Fig.  4.  (a)  Grid  tagging,  (b)  Fourier  spectrum  of  the  image  with 

the  corresponding  bandpass  filters,  (c)  Phase  images  obtained  from 
the  inverse  Fourier  Transform  of  the  bandpass  filtered  spectral  peaks 
added  together  (d)  Tag  locations  displayed  on  top  of  the  heart 


Fig.3.(c)  and  (d)  are  the  corresponding  Fourier  spectra 
of  these  images,  respectively.  A  bandpass  filter  is  applied 
to  the  deformed  tag  spectrum  to  extract  the  first  spectral 
peak.  The  IDFT  of  this  peak  is  a  complex  image.  The 
phase  map  of  this  complex  image  is  displayed  in  Fig.3.(e). 
The  agreement  between  the  tag  deformation  and  the  phase 
map  can  be  seen  clearly  in  this  image.  The  magnitude  of 
the  IDFT  of  this  peak  is  shown  in  Fig.3.(f).  As  expected 
from  Eq.7,  a  blur  in  the  magnitude  image  is  observed  be¬ 
cause  of  the  convolution  with  the  jinc  function.  Morpho¬ 
logical  closing  followed  by  a  simple  thresholding  is  applied 
to  the  image  in  Fig.3.(b)  to  segment  the  region  of  interest 
and  this  mask  is  applied  to  the  tag  detected  tag  points 
that  are  extracted  from  the  phase  image  using  a  simple 
thresholding  as  shown  in  Fig.3.(g).  Detected  tag  points 
are  added  on  top  of  the  image  in  Fig.3.(f)  to  demonstrate 
the  effectiveness  of  the  method  in  Fig.3.(h). 

B.  Grid  Tagging 

One  of  the  advantages  of  the  method  is  that  it  can  be 
easily  applied  to  two  dimensional  (grid)  tagging.  Separat¬ 
ing  the  grid  type  tags  into  two  orthogonal  components  us¬ 
ing  the  Fourier  Spectrum  is  first  proposed  by  Zhang  et  al. 
[10].  Mathematically,  rotation  in  spatial  domain  results  in 
an  equivalent  rotation  in  the  Fourier  domain.  Therefore, 
grid  tagging  is  simply  another  extra  step  of  90°  rotation 
plus  tagging  of  the  already  tagged  image.  This  will  create 
the  corresponding  spectral  peaks  in  orthogonal  directions. 

Fig. 4.  (a)  and  (b)  present  a  grid  tagged  image  and  its 
corresponding  Fourier  spectrum,  respectively.  The  user- 
defined  bandpass  filters  are  are  also  shown  for  the  first 
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spectral  peaks  in  each  tagging  direction.  Fig.4.(c)  shows 
the  resulting  phase  image  by  adding  the  two  phase  images 
obtained  from  the  IDFT  of  the  corresponding  bandpass 
filtered  spectral  peaks.  One  can  see  the  grid  pattern  in 
the  phase  image  indicating  the  agreement  of  phase  values 
and  the  tag  deformation.  From  these  two  phase  images  tag 
patterns  are  extracted  using  a  simple  thresholding  around 
neighborhood  of  n  radians.  Thresholded  images  are 
combined  with  a  logical  OR  operation.  Epicardial  con¬ 
tour  is  manually  traced  in  the  image  and  applied  as  a 
mask  to  remove  the  detected  tag  patterns  outside  the  re¬ 
gion  of  interest.  Resulting  image  is  added  on  top  of  the 
original  image  as  shown  in  Fig.4.(d)  to  display  the  agree¬ 
ment  between  the  detected  and  actual  tag  locations. 

IV.  Discussion 

One  of  the  issues  in  the  method  is  the  proper  bandpass 
filter  selection.  Currently,  the  user  is  asked  to  define  the 
band-pass  filter  by  drawing  a  circular  or  polygonal  region 
of  interest  on  the  Fourier  spectrum  of  the  image.  Apply¬ 
ing  a  Gaussian  roll  off  as  done  in  [1]  and  [7]  would  re¬ 
duce  the  ringing  artifact.  Another  improvement  would  be 
automating  the  filter  selection  process.  A  simple  thresh¬ 
olding  in  the  Fourier  Spectrum  image  would  give  us  the 
peak  locations  and  then  the  bandpass  filter  can  be  easily 
defined  without  user  intervention  according  to  Property  1. 
Currently,  the  same  filter  is  used  for  the  whole  image  set. 
However  the  DC  component  of  the  Fourier  spectrum  gets 
larger  in  time  because  of  the  T 1  recovery  of  the  underly¬ 
ing  tagged  tissue  and  might  interfere  with  the  band-pass 
filtered  peak.  An  adaptive  scheme  can  be  easily  imple¬ 
mented  to  update  the  extend  of  the  bandpass  filter  to 
eliminate  this  potential  problem.  However  the  drawback 
of  this  approach  is  that  reducing  the  extend  of  the  region 
also  means  limiting  the  ability  to  detect  deformed  tags  as 
deformations  result  in  deviations  from  the  center  of  the 
band-pass  filter  region. 

Tag  tracking  is  completed  under  a  minute  right  now 
including  the  bandpass  filter  selection  for  a  set  of  30  im¬ 
ages  with  256  x  256  dimensions  on  a  Ultra  Sparc  400  MHz 
processor. 


The  run  time  would  reduce  five  to  ten  seconds  once  the 
user  intervention  at  the  filter  selection  is  eliminated.  The 
algorithm  can  further  be  speeded  up  by  parallelizing  the 
process  as  the  tag  tracking  is  performed  independently  for 
each  frame. 

Our  emphasis  is  now  on  developing  a  fully  automated 
tag  analysis  program  using  this  phase  information  as  the 
tag  tracking  module. The  goal  is  to  complete  the  myocar¬ 
dial  contour  segmentation  and  data  analysis  under  4  min¬ 
utes  to  keep  the  whole  process  under  5  minutes  to  make 
it  useful  in  clinical  setting. 
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